function [fd, X, w] = tempfilter(d, opts)
% tempfilter  - filter temporal data
%
% FORMAT:       fd = tempfilter(d, opts)
%
% Input fields:
%
%       d           data to filter
%       opts        mandatory struct but with optional fields
%        .spat      enable spatial filtering (default: false)
%        .spkern    smoothing kernel in sampling units (default: [2, 2, 2])
%        .tdim      temporal filter dimension (default: 1)
%        .temp      enable temporal filtering (default: true)
%        .tempdct   DCT-based filtering (min. wavelength, default: Inf)
%        .tempdt    detrend (default: true, is overriden by dct/sc)
%        .templp    temporal lowpass (smoothing) kernel in units (def: 0)
%        .tempsc    sin/cos set of frequencies (number of pairs, def: 0)
%        .trobust   perform temporal filtering robustly (default: false)
%
% Output fields:
%
%       fd          filtered data (in input datatype, scaled if needed)
%       X           filtering matrix
%       w           weights of robust regression

% Version:  v0.7g
% Build:    9031016
% Date:     Mar-10 2009, 4:33 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 2 || ...
   (~isnumeric(d) && ...
    ~istransio(d)) || ...
   ~isstruct(opts) || ...
    numel(opts) ~= 1
    error( ...
        'BVQXtools:BadArgument', ...
        'Invalid or missing argument in call.' ...
    );
end
if ~isfield(opts, 'tdim') || ...
   ~isa(opts.tdim, 'double') || ...
    numel(opts.tdim) ~= 1 || ...
   ~any(opts.tdim == (1:ndims(d)))
    opts.tdim = 1;
end
if numel(d) == size(d, opts.tdim)
    opts.spat = false;
end
if ~isfield(opts, 'spat') || ...
   (~islogical(opts.spat) && ...
    ~isnumeric(opts.spat)) || ...
    numel(opts.spat) ~= 1
    opts.spat = false;
else
    opts.spat = (opts.spat ~= 0);
end
if ~isfield(opts, 'temp') || ...
   (~islogical(opts.temp) && ...
    ~isnumeric(opts.temp)) || ...
    numel(opts.temp) ~= 1
    opts.temp = true;
else
    opts.temp = (opts.temp ~= 0);
end
if opts.spat
    if ~isfield(opts, 'spkern') || ...
       ~isa(opts.spkern, 'double') || ...
        numel(opts.spkern) ~= 3 || ...
        any(isinf(opts.spkern) | isnan(opts.spkern) | opts.spkern < 0 | opts.spkern > 10)
        opts.spkern = [2, 2, 2];
    else
        opts.spkern = opts.spkern(:)';
    end
    if all(opts.spkern) < 0.5
        opts.spat = false;
    end
end
if opts.temp
    if ~isfield(opts, 'tempdct') || ...
       ~isa(opts.tempdct, 'double') || ...
        numel(opts.tempdct) ~= 1 || ...
        isinf(opts.tempdct) || ...
        isnan(opts.tempdct) || ...
        opts.tempdct < 4
        opts.tempdct = Inf;
    end
    if ~isfield(opts, 'tempdt') || ...
       ~islogical(opts.tempdt) || ...
        numel(opts.tempdt) ~= 1
        opts.tempdt = true;
    end
    if ~isfield(opts, 'templp') || ...
       ~isa(opts.templp, 'double') || ...
        numel(opts.templp) ~= 1 || ...
        isinf(opts.templp) || ...
        isnan(opts.templp) || ...
        opts.templp < 0 || ...
        opts.templp > 12
        opts.templp = 0;
    end
    if ~isfield(opts, 'tempsc') || ...
       ~isa(opts.tempsc, 'double') || ...
        numel(opts.tempsc) ~= 1 || ...
       ~any(opts.tempsc == (1:floor(0.5 * (size(d, opts.tdim)-1))))
        opts.tempsc = 0;
    end
    if ~isfield(opts, 'trobust') || ...
       ~islogical(opts.trobust) || ...
        numel(opts.trobust) ~= 1
        opts.trobust = false;
    end
    if isinf(opts.tempdct) && ...
       ~opts.tempdt && ...
        opts.templp == 0 && ...
        opts.tempsc == 0
        opts.temp = false;
    elseif opts.tempsc > 0
        opts.tempdct = Inf;
        opts.tempdt = false;
    end
end

% return if yet nothing to do...
X = zeros(size(d, opts.tdim), 0);
if ~opts.spat && ...
   ~opts.temp
    fd = d;
    return;
end

% get data
dt = class(d);
if istransio(d)
    if ~strcmp(dt, 'double')
        fd = single(resolve(d));
    else
        fd = resolve(d);
    end
elseif ~strcmp(dt, 'double')
    fd = single(d);
else
    fd = d;
end

% permute if necessary
di = opts.tdim;
if di > 1
    tperm = 1:ndims(fd);
    tperm = [di, tperm(tperm ~= di)];
    [sperm, fperm] = sort(tperm(:));
    fd = permute(fd, tperm);
end

% get original sizing and number of time points
vs = size(fd);
nv = vs(1);

% DCT overrides detrending
if opts.temp && ...
   ~isinf(opts.tempdct)
    opts.tempdct = floor(2 * nv / opts.tempdct);
    opts.tempdt = false;
end

% spatial filter first if needed
if opts.spat
    
    % if too few dims, make it so
    fds = size(fd);
    if numel(fds) < 4
        fd = reshape(fd, [fds(1), ones(1, 4 - numel(fds)), fds(2:end)]);
        opts.spkern(1:(4 - numel(fds))) = 0;
    end
    ssr = size(fd);
    ssf = ssr(2:end);
    ssr(1) = 1;
    
    % use smoothdata3 to iterate over volumes
    for vc = 1:nv
        fd(vc, :, :, :, :) = reshape( ...
            smoothdata3(reshape(fd(vc, :, :, :, :), ssf), opts.spkern), ssr);
    end
    
    % reshaping again?
    if numel(fds) < 4
        fd = reshape(fd, fds);
    end
end

% temporal filter
if opts.temp
    
    % reshape if necessary
    if numel(vs) > 2
        fd = reshape(fd, nv, prod(vs(2:end)));
    end

    % only detrend
    if opts.tempdt
        
        % build X
        X = [ones(nv, 1), (-1:(2/(nv-1)):1)'];
        
    % DCT-based filtering
    elseif ~isinf(opts.tempdct)
        
        % build X
        X = [ones(nv, 1), zeros(nv, opts.tempdct)];
        n = 0:(nv - 1);
        for dc = 1:opts.tempdct
            X(:, dc + 1) = cos(pi * (2 * n + 1) * dc / (2 * nv));
        end
        
    % sin/cos set filtering
    elseif opts.tempsc > 0
        
        % build X
        X = [ones(nv, 1), zeros(nv, opts.tempsc)];
        n = 0:(nv - 1);
        for pc = 1:opts.tempsc
            X(:, (pc * 2)) = sin(pi * (2 * n + 1) * pc / nv);
            X(:, (pc * 2 + 1)) = cos(pi * (2 * n + 1) * pc / nv);
        end
    end
    
    % perform regression
    if opts.trobust
        if nargout < 3
            b = fitrobustbisquare_img(X, fd')';
        else
            [b, w] = fitrobustbisquare_img(X, fd');
            b = b';
            w = reshape(w', vs);
            if di > 1
                w = permute(w, fperm');
            end

        end
    else
        b = pinv(X' * X) * X' * double(fd);
        if nargout > 2
            w = ones(vs);
        end
    end
    
    % keep mean
    b(1, :) = 0;
    
    % subtract from data
    fd = fd - X * b;
    
    % low-pass filter?
    if opts.templp > 0
        
        % get kernel
        k = smoothkern(opts.templp, 1e-5);
        sk = zeros(size(fd));
        ks = zeros(nv, 1);
        khs = 0.5 * (numel(k) - 1);
        
        % sum up
        for kc = -khs:khs
            sf = max(1, 1-kc);
            st = min(nv, nv-kc);
            tf = max(1, 1+kc);
            tt = min(nv, nv+kc);
            
            sk(tf:tt, :) = sk(tf:tt, :) + k(1 + khs + kc) * fd(sf:st, :);
            ks(tf:tt) = ks(tf:tt) + k(1 + khs + kc);
        end
        
        % rebuild fd
        fd = sk;
        sk = [];
        ks = 1 ./ ks;
        fd = fd .* ks(:, ones(1, size(fd , 2)));
    end
    
    % reshape again
    fd = reshape(fd, vs);
    
    % remove mean from filter matrix
    X(:, 1) = [];
end

% datatype
if ~strcmp(dt, class(fd))
    
    % get minmaxmean
    mmm = minmaxmean(fd);
    rng = eps + (mmm(2) - mmm(1));
    switch(dt)
        case {'int16'}
            if mmm(1) < -32768 || ...
                mmm(2) >= 32767.5
                fd = -32768 + (65535 / rng) * (fd - mmm(1));
            end
        case {'int32'}
            if mmm(1) < -2147483648 || ...
                mmm(2) >= 2147483647.5
                fd = -2147483648 + (4294967295 / rng) * (fd - mmm(1));
            end
        case {'int8'}
            if mmm(1) < -128 || ...
                mmm(2) >= 127.5
                fd = -128 + (255 / rng) * (fd - mmm(1));
            end
        case {'uint16'}
            if mmm(1) < 0 || ...
                mmm(2) >= 65535.5
                fd = (65535 / rng) * (fd - mmm(1));
            end
        case {'uint32'}
            if mmm(1) < 0 || ...
                mmm(2) >= 4294967295.5
                fd = (4294967295 / rng) * (fd - mmm(1));
            end
        case {'uint8'}
            if mmm(1) < 0 || ...
                mmm(2) >= 255.5
                fd = (255 / rng) * (fd - mmm(1));
            end
    end
    eval(['fd=' dt '(round(fd));']);
end

% re-permute if necessary
if di > 1
    fd = permute(fd, fperm');
end
